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Protected areas are a key tool in the conservation of global biodiversity and carbon stores. We conducted a global test of the 
degree to which more than 18,000 terrestrial protected areas (totalling 5,293,217 km?) reduce deforestation in relation to 
unprotected areas. We also derived indices that quantify how well countries’ forests are protected, both in terms of forested 
area protected and effectiveness of protected areas at reducing deforestation, in relation to vertebrate species richness, 
aboveground forest carbon biomass and background deforestation rates. Overall, protected areas did not eliminate deforesta- 
tion, but reduced deforestation rates by 41%. Protected area deforestation rates were lowest in small reserves with low back- 
ground deforestation rates. Critically, we found that after adjusting for effectiveness, only 6.5%—rather than 15.7%—of the 
world's forests are protected, well below the Aichi Convention on Biological Diversity’s 2020 Target of 17%. We propose that 
global targets for protected areas should include quantitative goals for effectiveness in addition to spatial extent. 


pre-human background rates, suggesting that we are in the 
midst of a sixth mass extinction event'. Diverse efforts to 
conserve biodiversity include captive breeding programmes, legal 
protections for individual species, restrictions on wildlife trade, 
control of invasive species and the establishment of protected areas” 
(PAs). PAs are a cornerstone of many conservation programmes 
and can be effective in reducing overexploitation, habitat loss and 
many other threats within their boundaries’. The percentage of 
Earth's terrestrial area that is protected has increased substantially 
over the past decade and is now approaching the Aichi Convention 
on Biological Diversity’s 2020 Target of 17%*°. However, PA cov- 
erage alone is an inadequate conservation metric because nearly 
one-third of protected land is under intense human pressure‘ and 
‘paper parks’ prove ineffective’. Moreover, prioritizing area pro- 
tected over other metrics can indirectly lead to declines in effec- 
tiveness’. Assessing the effectiveness of PAs is challenging because 
biodiversity responses to protection are often difficult to measure’. 
One of the most basic functions of PAs is to provide habitat for 
species. Most of Earth's terrestrial species rely on forest habitat’. 
Consequently, we focus here on deforestation as a form of habi- 
tat loss, although other types of habitat loss can be important in 
non-forest ecosystems. Deforestation increases the extinction risk of 
forest-associated species’ and reduces the amount of carbon seques- 
tered into organic biomass’. Thus, the extent to which PAs limit 
deforestation is often a key component of their effectiveness, both 
in terms of species conservation and carbon storage services. Unlike 
species populations”'»”’, forest cover and change can be consistently 
mapped globally at high resolution using remote sensing, which 
makes such data well-suited to global analyses'*. Previous research 
has demonstrated that PAs conserve forest habitat’, that mixed-use 
and indigenous PAs can be more effective than strictly protected 
PAs'*" and that PAs were most effective in preventing deforesta- 
tion in Australasia and least effective in Asia, with effectiveness 


( urrent species extinction rates are ~1,000 times higher than 


positively associated with countries’ gross domestic product (GDP) 
per capita’®. This parallels the finding that PA effectiveness at limit- 
ing human pressure increases was positively associated with human 
development index”. While PA establishment can lead to increased 
deforestation rates in nearby areas'* (‘leakage’), the opposite out- 
come (‘blockage’) is much more common globally”. 

Here, we build on this previous research to conduct the first (to 
the best of our knowledge) comprehensive, global analysis of the 
effectiveness of PAs with respect to limiting forest loss. We defined 
PA effectiveness (with regard to limiting deforestation) based on 
deforestation rates within PAs compared with rates in matched con- 
trol areas with similar characteristics. We modelled PA deforesta- 
tion rates while controlling for background rates in matched control 
areas using a diverse set of predictors, including: nearby population 
densities; reserve size, age and management category; and GDP per 
capita. As geographical variation and local or regional context can 
affect how these predictors relate to deforestation rates, we adopted 
a spatially and non-spatially varying coefficient (SNVC) model- 
ling approach, which relaxes the assumption of constant effect 
sizes. By allowing relationships to vary geographically, our analysis 
can shed light on differences in results among previous PA 
effectiveness studies. 

From a conservation perspective, regions with high species rich- 
ness and carbon storage should ideally have high PA cover and 
effectiveness, because protecting biodiversity and carbon stocks can 
be important functions of PAs”. In addition to our primary model- 
ling effort, we derived a national scale index of effective area pro- 
tected, which we compared with total forest vertebrate richness and 
carbon stocks across all countries on the global scale. This allows 
for the quantification of how well effective protection aligns with 
biodiversity and carbon sequestration”', and provides a transpar- 
ent assessment of which countries are the most under-protected 
given their biodiversity and carbon stocks, and which countries 
have excelled at effective forest protection. This comparative, global 
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Fig. 1 | Locations of the 18,171 PAs in our main analysis and global deforestation rates. a, PAs are grouped based on their IUCN categories (cat.): Strict (I- 
IV), Nonstrict (V-VI) and Unknown. Most of the PAs are located within regions of higher GDP per capita, especially the eastern United States and Europe. 
However, many of these PAs are very small. b, Estimated annual deforestation (tree cover loss) rates are for the period 2001 to 2018. Including PAs from 
all of Earth's forested regions enabled our analysis to assess geographic variation in effectiveness on the global scale. 


assessment complements the current coverage target of protection 
associated with the Aichi Biodiversity Target 11 and other global 
effectiveness evaluations. 

While the subject of PA effectiveness has received considerable 
attention in the literature, most previous ecological evaluations of 
PAs have been regional in scope, which precludes testing global 
hypotheses about political and economic predictors of PA effective- 
ness. However, there are examples of global studies using either the 
change in human pressure®”” or biodiversity metrics”. Our use 
of a global set of PAs and an SNVC modelling framework allowed 
us to test the relative performance of all forested PAs across the 
full gradients of latitude, population density, GDP per capita and 
other important predictor variables, while accounting for regional 
scale variability. Furthermore, our SNVC modelling approach 
permits testing of whether the effects of these variables vary geo- 
graphically or can be treated as constant. Finally, we provide what 
we believe is the first assessment of the degree to which PA effec- 
tiveness and coverage are congruent with areas in greatest need 
of protection, including those with the highest biodiversity and 
carbon stocks. We thus provide the first (to the best of our knowl- 


edge) global quantitative estimates of the most under-protected 
countries on Earth. 


Results 

Prior to matching with control areas, our primary PA data- 
set contained 25,348 PAs. However, 7,177 (28.3%) of these PAs 
could not be matched with any unprotected pixels having similar 
matching-covariate values. Consequently, after applying coarsened 
exact matching to identify control areas, our final dataset contained 
18,171 PAs, with a total area of 5,293,217 km? (Fig. 1 and Extended 
Data Fig. 1). Overall, PAs reduced, but did not eliminate, deforesta- 
tion; the median annual deforestation rate in control areas (0.54%; 
s.d.=2.21%) was 4.97 times higher than within PAs (0.11% per 
year; s.d. =2.45%) (Fig. 2). In absolute terms, the 18,171 PAs in our 
analysis had an average annual forest loss rate of ~1.53 Mha. In our 
analysis, 28.7% of the PAs did not have any forest loss. Among PAs 
with known management category, deforestation rates were high- 
est in nonstrict PAs in Africa (0.31% per year), Europe (0.29% per 
year) and South America (0.19% per year), and lowest in strict PAs 
in Oceania (0.02% per year) (Fig. 3). General patterns in PA forest 
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Fig. 2 | Deforestation rate distributions. Histograms of deforestation rates within PAs (left column) and associated control areas (right column) grouped 
by IUCN category (rows) for the 18,171 PAs in our main analysis. Deforestation rates (2001-2018) are expressed as mean annual rate relative to tree 
cover in 2000. Deforestation rates tended to be greater in control areas, providing evidence of PA effectiveness. Deforestation rate was truncated at 5% 


(covering more than 98% of the data) and log(1+ x) transformed for plotting. 


loss for the different continents and levels of protection tended to 
be similar when net forest loss was used instead of total forest loss 
(Fig. 3 and Extended Data Fig. 2). 

We identified 9,875PAs established between 2002 and 2017 
that were suitable for inclusion in our spatiotemporal analysis. 
The establishment of these PAs was associated with a moderate 
increase in the deforestation rate (average = 0.19%; s.e.m.=0.02%), 
whereas control areas saw a larger increase in the deforestation 
rate (average=0.61%; s.e.m.=0.02%) over the same time span 
(Extended Data Fig. 3). This overall pattern was observed for 
both strict and nonstrict PAs in lower and higher GDP countries 
(Extended Data Fig. 3). 


Predictors of PA effectiveness. We found all predictors of PA defor- 
estation rates to be spatially varying with the exceptions of travel 
time to nearest densely populated area (estimate = 0.066; standard 
error=0.071; P=0.356), threatened forest species richness (esti- 
mate = 0.080; standard error =0.037; P= 0.029) and non-threatened 
forest species richness (estimate=0.018; standard error = 0.004; 
P<0.001) (Fig. 4 and Extended Data Figs. 4 and 5). Reserve area 
and background (that is, control area) deforestation rate were both 
generally positively associated with PA deforestation, although, in 
both cases, the effects appeared to be stronger in high latitudes and 
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weaker in the tropics (Fig. 4). The effects of the other predictors 
were often relatively weak or inconsistent (Extended Data Fig. 4). 


National level PA effectiveness scoring. Overall, the mean defores- 
tation rate within PAs was 41.1% lower than in control areas (0.62% 
per year compared with 1.05% per year). This is analogous to defor- 
estation within PAs occurring at background rates over 58.9% of 
their area and at a rate of 0% over 41.1% of their area, on average 
(0.0062 ~0.411X0+0.589 x0.0105). Thus, the 17% Aichi 
Target goal (with 0% PA forest loss) is equivalent to 41.3% of 
land protected after accounting for deforestation within PAs (as 
0.413 x 0.411 + 0.17)—a more than twofold increase. 

Globally, 15.7% of forest is formally protected; however, after 
adjusting for deforestation within reserves, this was reduced to only 
6.5%. That is, the 15.7% of forest protected with current defores- 
tation rates would have the same total deforestation rate as 6.5% 
protected with no deforestation and 9.2% protected with the con- 
trol area mean deforestation rate (0.157 X 0.0062 ~ 0.065 x0 + 0.09 
2x0.0105). Among the 63 countries that we considered, 34 (54%) 
have at least 17% of their forested area protected (Supplementary 
Table 1). However, countries varied greatly in terms of area pro- 
tected after adjusting for effectiveness, where effectiveness is 
defined as the ratio of control area deforestation to PA deforestation 
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Fig. 3 | Forest loss in and around PAs. Results are grouped by geographic region and PA IUCN category. For each PA, the variables shown are PA and 
control area forest loss. Forest loss is for the period 2001-2018 and is expressed as the annual deforestation rate relative to forest cover in 2000. Points 
correspond to median (across PAs) percentage forest loss. Error bar end points are the 1st and 3rd quartiles for this variable. Forest loss within PAs has 
generally been less than in nearby unprotected areas. The overall pattern of forest loss being lower in PAs than in control areas was generally consistent 


across both PA type and geographic region. 


(Fig. 5 and Supplementary Table 1). These adjusted area protected 
values ranged from 0.03 (China) to 2.36 (New Zealand) (Fig. 5). 
South Africa had the highest effectiveness score, 8.10, indicating 
that PAs there collectively had a greater than eightfold decrease in 
deforestation rates relative to matched control areas (Supplementary 
Table 1). Strikingly, there were no countries with both very high 
PA coverage (adjusted for effectiveness) and high forest species 
richness (Fig. 5). 

There was similarly high variation in the species threat index 
(ratio of forest vertebrate species richness to level of protec- 
tion after accounting for effectiveness), with values for the seven 
countries with at least 500 forest vertebrate species ranging from 
1,516 (Mexico) to 6,913 (Indonesia). The continents with the high- 
est mean adjusted threat indices were South America (2,571), 
Asia (2,196) and Oceania (1,834). The countries with the high- 
est deforestation rates were Sierra Leone (2.40% per year; 5.63% 
of forested area protected), Malaysia (2.16% per year; 6.29% pro- 
tected) and Cambodia (1.90% per year; 25.85% protected) (Fig. 6 
and Supplementary Table 1). Of the three countries with the great- 
est amounts of aboveground forest carbon, Russia had the highest 
threat index for carbon (109), followed by the United States (57) and 
Brazil (29) (Fig. 6 and Supplementary Table 1). The 11 countries 


in our analysis with less than 10% of their forested area protected 
had a total of 72.6 Gt aboveground forest carbon (Supplementary 
Table 1). Forest carbon and total forest species richness were 
moderately correlated (r=0.37). 


Diagnostic and sensitivity analyses. Overall quality of matching 
was high, with absolute standardized biases ranging from 0.004 
for slope to 0.038 for population density. The Rosenbaum bounds 
excluded zero up to '=4, suggesting that our main result about PA 
effectiveness relative to control areas is at least moderately robust to 
hidden biases. When stricter matching criteria were applied (9-10 
classes per continuous matching covariate), the resulting dataset 
size decreased from 18,171 observations to 13,291 observations. For 
this new dataset, the mean deforestation rate in PAs was estimated 
to be 42.7% lower than in control areas (compared with 41.1% with 
our main dataset). Deforestation rate patterns by level of protection 
and geographic region were generally similar (Extended Data Fig. 
6). Overall modelling conclusions about the effects of background 
rate deforestation and reserve size were also similar, although the 
effects of strict protection, reserve age and population density 
(along with travel time) were found to not be spatially varying for 
the smaller dataset (Extended Data Fig. 7). 
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Fig. 4 | Effects of control area deforestation rate and PA area on PA deforestation rates. The maps show spatially varying coefficient model results 
(coefficient estimates, standard errors and false-discovery-rate-adjusted P values) for these predictors of PA deforestation. These results indicate relatively 
consistent positive relationships for control area (background) deforestation and reserve area. Extended Data Fig. 4 shows results for the complete set of 
spatially varying predictors. Only coefficients with associated P values less than 0.05 are mapped. 


Discussion 

PAs have been put forth as an important policy tool to ensure global 
biodiversity conservation and carbon storage in the face of expand- 
ing human resource exploitation’. However, limitations to enforce- 
ment and monitoring can reduce their effectiveness”. Indeed, 
our results indicate that PAs are rarely, if ever, strictly ‘protected 
from deforestation, but rather have the effect of slowing forest loss 
in relation to matched control areas (Fig. 3). This is consistent with 
underfunding being common, especially in countries with lower 
GDP per capita where many reserves may lack the equipment, staff 
and resources needed for effective management***’. Importantly, 
while prevention of deforestation is just one of many potential PA 
management goals, it may be correlated with other metrics of PA 
management success. For example, forest species are expected to 
have reduced likelihood of population declines in PAs with lower 
deforestation rates, either owing to greater forested habitat availabil- 
ity, or other benefits of more effective management, such as effec- 
tive suppression of illegal hunting. The correlation, albeit moderate, 
between forest biodiversity and carbon stocks provides further 
evidence of potential co-benefits, but also suggests there may be 
important tradeoffs that must be considered when selecting loca- 
tions for reserves (Supplementary Table 1). 

Earth's forests are a major carbon sink, with an estimated uptake of 
2.3 Pg Cyr“ between 2000 and 2007”. For comparison, the estimated 
total emissions rate over this period was 8.7PgCyr', of which 
roughly 1.1 PgCyr7! was due to tropical land use change”. As PAs 
have the potential to limit deforestation, they can be of use in miti- 
gating future climate change as well as protecting biodiversity”. 
Our finding that the effects of PAs on forest loss are highly variable 
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has substantial implications for the implementation of carbon pay- 
ment systems such as the United Nations’ Reducing Emissions from 
Deforestation and Forest Degradation (REDD+) programme and 
for future payment for environmental services schemes in that there 
is a tradeoff between focusing funding on areas where PAs are com- 
mon and effective (but may have less room for improvement) and 
areas where PAs are less common and less effective”. Although the 
actual implementation and impact of REDD+ has been hampered 
bya lack of commitment from the potential funding sources in coun- 
tries with higher GDP per capita™, several countries are construct- 
ing specific payment for environmental services schemes supported 
by national funds’. Such funds could make use of regional threat 
indices similar to those in our analysis (Supplementary Discussion), 
allowing them to identify the most vulnerable hotspots within their 
national PAs and prioritize actions in parks that will have the most 
effective outcomes. Using these fine-scale indices, internal funding 
could be directed towards areas with abundant biodiversity, high 
forest carbon stocks and limited protection. 

Larger PAs tended to have higher forest loss rates, even after 
accounting for loss in paired control areas, although this difference 
was lower in tropical regions (Extended Data Fig. 4). This is con- 
sistent with previous work that has shown that both budgets and 
staffing on a per-area basis are negatively correlated with reserve 
size”, although costs also have been shown to decline rapidly with 
increasing reserve size. Thus, PA size may be an especially impor- 
tant consideration for biodiverse countries with limited protection, 
which could benefit from both greater area protected and more 
effective protection (Fig. 5 and Supplementary Table 1). Our finding 
on PA size may be partly attributable to controlling for background 
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Fig. 5 | Forest biodiversity threat index. a, Number of forest vertebrate species versus area protected adjusted (that is, multiplied) by country-level PA 
effectiveness (based on forest loss in PAs compared with forest loss in matched control areas). The colours indicate the ratio of these variables, which we 
term species threat index, and provide insight into which countries have exceptionally high forest species richness relative to their level of protection. 

b, Threat index map showing only countries with at least 15 PAs in our main analysis, at least 5 forest obligate vertebrates and at least 10,000 km? forest. 
The adjusted threat index tended to be greatest in the tropics as expected, given the large number of species in these regions. 


deforestation rates given that large reserves may tend to be in areas levels of use allowed within nonstrict PAs, particularly when they 
with relatively little human impact and thus lower deforestation overlap communal lands. 

pressure. While strict protection was more effective than nonstrict After adjusting for effectiveness at limiting deforestation, 
protection in a few small regions, the effect was not spatially robust New Zealand scored highest for forest protection (Fig. 5 and 
(Extended Data Fig. 4). This may be a consequence of differences in Supplementary Table 1). In this country, widespread loss of 
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Fig. 6 | Annual deforestation rate and total aboveground forest carbon versus adjusted forested area protected. a,b, Annual deforestation rate (a) and 
total aboveground forest carbon (b) versus adjusted forested area protected. The adjustment was made by multiplying the proportion of forested area 
protected by country-level median PA effectiveness (based on forest loss near PAs compared with forest loss within PAs). The colours indicate the ratios 
of these variables, which we term threat indices (for forest loss and carbon), and provide insight into which countries have exceptionally high deforestation 
rates or forest carbon stocks relative to their levels of protection. The dataset is provided in Supplementary Table 1. DR Congo, Democratic Republic of 


the Congo. 


indigenous land cover occurred in recent years, but most of this 
loss (93.9%) occurred outside PAs”. This suggests that PAs in New 
Zealand may be particularly effective at preventing deforestation. 
New Zealand also had the second highest PA effectiveness score 
(7.53), after South Africa (8.10). However, these results provide only 
a partial picture because major forest loss occurred in New Zealand 
prior to 2000”. Conversely, among countries with at least 1,000 for- 
est vertebrate species, Indonesia is notable in that it had both the 
least amount of forested land protected and the highest species 
threat index score (Fig. 5). Biodiversity in Indonesia is under threat 
owing to this country’s high rate of deforestation, increasing oil 
palm production and illegal wildlife trade“. Furthermore, Southeast 
Asia has the greatest proportions of endemic bird and mammal 
species, and the highest rate of forest loss”. 


Conservation implications. The species-based threat index that 
we present quantifies countries’ forest vertebrate biodiversity rela- 
tive to their levels of protection (Fig. 5 and Supplementary Table 
1). It provides what we believe is the first quantitative measure of 
the extent to which forest biodiversity and protection align. The 
finding that many countries with higher GDP per capita are com- 
paratively well protected given their levels of forest biodiversity 
(Fig. 5) is not surprising considering the financial costs of PAs, and 
suggests that economic growth could eventually allow for greater 
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protection. Unfortunately, environmental degradation and species 
loss often occur in the early periods of development“, and it may 
be more difficult to reintroduce native species than to increase PA 
area or effectiveness. Consequently, it should not be assumed that 
countries with lower GDP per capita can readily restore biodiversity 
in the long term after income levels increase and forest protection 
measures are enhanced. Similarly, existing forests store substantial 
amounts of carbon (Fig. 6), and their loss may not be easily com- 
pensated for (for example, through reforestation), especially over 
short timescales. These timescale differences demonstrate the 
urgent need to expand and strengthen PA networks where they can 
be most beneficial. 

The Convention on Biological Diversity’s global PA target (Aichi 
Target 11) sets a global, rather than country-specific, percentage 
goal (17%). While this goal often motivates countries targeting 17% 
protection, based on our threat index, different goals for individual 
countries could instead be advocated, with higher targets in coun- 
tries where biodiversity is most threatened*. To achieve these new 
goals, we therefore see it as vital that countries with higher GDP per 
capita provide support for the establishment and maintenance of 
PAs necessary to achieve these higher targets in biodiverse countries 
with lower GDP per capita, ensuring distributional equity in global 
conservation targets. It is critical that such support is provided 
in an equitable way that promotes social justice and sustainable 
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development, and is not used to justify ignoring the protection of 
less biodiverse areas globally”. While proximate threats to biodiver- 
sity tend to be local, the drivers of these threats are often global”. 
Similarly, many potential benefits and co-benefits associated with 
biodiversity are accrued on both local and global scales. Thus, the 
country-specific threat indices reported here highlight human- 
ity’s failure to allocate conservation funding, research and other 
resources where they are most needed. 

Many biodiverse countries with substantial forest carbon stocks 
are unable to achieve effective conservation in isolation—especially 
given outside demands for their natural resources. One way to 
reduce these demands is through the establishment of multinational 
import bans on deforestation-associated commodities“, which 
must be coupled with internal enforcement and more efficient use 
of forest resources (for example, through technological develop- 
ments). Additionally, strategies such as environmental certification 
and labelling schemes could be used to give consumers in wealthy 
countries an opportunity to reward environmentally friendly prac- 
tices“. This may be especially effective in reducing the occurrence 
of selective illegal logging, which is often difficult to monitor and 
may be prevalent in large reserves that tend to have high apparent 
deforestation rates and are more likely to be downsized, down- 
graded or degazetted* (Fig. 4). More research is needed to assess 
the possible benefits and drawbacks of import bans and certifica- 
tion and labelling schemes—a topic beyond the scope of our pres- 
ent analysis. Although global conservation problems require global 
solutions, local factors must also be considered, including how 
conservation programmes can best make use of local knowledge, 
benefit and empower communities, and contribute to long-term 
economic development and poverty reduction”. 

As of 2018, 14.9% of land is protected, which is close to the 17% 
Aichi goal“. However, no similar numeric target has been stated 
for PA quality. This creates a policy incentive to value total PA 
area above PA effectiveness”. Unfortunately, deforestation rates in 
large PAs (which cover a sizeable area) seem to be relatively high 
(Fig. 4). To date, more than 40 approaches have been developed to 
assess reserve context, management inputs and design”, involving 
techniques including questionnaires”, selection of focal conserva- 
tion targets”? and engagement with stakeholders”. However, such 
assessments do not directly capture outcomes. Thus, estimates of 
PA deforestation rates relative to those in matched control areas rep- 
resent a simple, outcome-focused complement to these methods, 
acting as a useful starting point for measuring effectiveness given 
the many co-benefits of reducing deforestation and growing for- 
ests to their ecological potential, including increased carbon stor- 
age, improved water and air quality, recreation opportunities, and 
reduced erosion®™. While we assessed PA deforestation rates rela- 
tive to matching control areas, it is important to note that absolute 
deforestation rates can also be of interest. This is especially true in 
cases where background deforestation rates are exceptionally high, 
potentially leading to PAs with moderate deforestation rates being 
identified as highly effective relative to control areas. 

By adjusting area protected using a measure of PA effectiveness, 
we have combined PA quantity (area) and quality (effectiveness) 
into a single metric (Figs. 5 and 6). This combination represents a 
small step towards preventing ‘perverse outcomes, wherein man- 
agers optimize PAs using a single, imperfect metric such as total 
area protected’, and can complement existing PA effectiveness 
assessments”'”**. Taken together, the combination of PA area and 
PA effectiveness can be used to set more robust targets for PAs—an 
important consideration for the upcoming Fifteenth meeting of the 
Conference of the Parties to the Convention on Biological Diversity 
conference, where new biodiversity conservation targets will be set. 
Scientists have called for half of Earth’s surface to be protected”. 
Although such direct statements can attract widespread public sup- 
port, the success or failure of future protection targets to conserve 


biodiversity and carbon stores will depend on whether or not they 
include sub-targets or goals that can be readily measured, connected 
to species and populations, and achieved in a way that benefits both 
humans and the environment on which we depend. 


Methods 


We estimated forest change rates within terrestrial PAs using the World Database 
on Protected Areas” and a set of global forest change maps”. We omitted PAs for 
which forest loss could not be accurately estimated, including recently established 
PAs, exceptionally small PAs and PAs outside forest biomes, using a series of 
filtering steps (Supplementary Methods). We resampled the tree cover (year 
2000), forest loss (2001-2018) and forest gain (2000-2012) maps from ref. ° to 
1km xX 1km resolution (Supplementary Methods). We then used forest loss to 
calculate an approximate annual deforestation rate and used forest loss and gain 
together to estimate net deforestation (Supplementary Methods). 

We identified associated control areas for PAs using a variant of 1-k coarsened 
exact matching” and eight ‘matching covariates’: elevation, slope, tree cover, 
travel time to nearest densely populated area™, population density”, country, 
ecoregion® and primary driver of forest cover loss. Our matched-comparison 
dataset consisted of one observation for each PA (treatment observations) 
and one observation for each 1-km raster pixel more than 10 km from all PAs 
(control observations). We avoided selecting control areas near PAs because 
these regions may have elevated deforestation rates due to local scale leakage 
(that is, displacement of deforestation from PAs to nearby areas as a result of 
protection)'**” and can differ substantially from adjacent PAs, making them 
unsuitable for use as counterfactual scenarios”. We first coarsened the combined 
dataset by discretizing continuous variables (Supplementary Methods) and then 
calculated the associated L1 imbalance indicating the extent to which treatment 
and control observations differ with respect to the matching covariates”. We then 
coarsened the dataset again and paired each PA with the set of control observations 
that had the same coarsened covariate values. Finally, we estimated forest loss rates 
within the control areas. 

To verify that similar treatment and control units were successful matched, we 
conducted diagnostic and sensitivity analyses”. We first calculated the standardized 


bias for each covariate, which is defined as (X, — X-)/o;, where X; is the treatment 


mean, X, is the control mean and g, is the treatment standard deviation“. Absolute 
standardized biases greater than 0.25 indicate that the quality of matching may 
be poor“. While standardized bias quantifies similarity with respect to observed 
covariates, there can also be hidden biases owing to unobserved confounding 
variables”. We explored sensitivity with respect to potential hidden biases using 
Rosenbaum bounds that quantify how results vary with respect to the odds of 
assignment (treatment or control) depending on unobserved variables”. We did 
this by calculating Rosenbaum bounds for the Hodges~Lehmann point estimate 
using the ‘hlsens’ function in the ‘rbounds’ R package” where deforestation rates 
were log(x+ 10~’) transformed to parallel subsequent modelling analyses. The odds 
of differential assignment due to variables not used for matching, I’, was allowed to 
vary from 1 to 6 in increments of 1. As a final sensitivity analysis, we recalculated 
our main results after matching using stricter criteria (Supplementary Methods). 

As additional modelling covariates, we determined threatened and 
non-threatened forest vertebrate species richness using International Union for 
Conservation of Nature (IUCN) Red List species range maps (Supplementary 
Methods). We also obtained year of establishment (where available), area and 
associated GDP per capita (country level; purchasing power parity, current 
international $) for each PA”. For each country, we used the average GDP per 
capita between 2000 and 2018 (excluding years without data). 

The most rigorous test of PA effectiveness is to compare deforestation 
between pre- and post-PA establishment in relation to associated control areas. 
Effectively, this is a ‘before-after, control-impact’ design’’. Unfortunately, most 
PAs were established before the origin of global deforestation data, precluding the 
use of this framework for our primary analysis. Nevertheless, nearly 10,000 PAs 
meeting our criteria for inclusion were established between 2002 and 2017, which 
enabled us to adopt this rigorous design for a subset of PAs as a secondary analysis 
(Supplementary Methods). 


Statistical modelling and hypotheses. Because the effects of predictors of 

PA deforestation may vary spatially due to differences in regional context, we 
adopted an SNVC modelling approach’”~” (Supplementary Methods). We used 
‘PA deforestation rate’ as the response variable and statistically controlled for 
deforestation rates in control areas by including this variable as a covariate. We 
included six other predictors of deforestation rates: population density, travel 

time to nearest densely populated area, PA age, GDP per capita, management 
category (Strict: IUCN category I-IV; Nonstrict: IUCN category V-VI) and PA 
area (Supplementary Methods). We formed the following set of a priori hypotheses 
based on the PA literature: 


1. Strictly protected PAs have lower deforestation rates than nonstrictly pro- 
tected PAs because they are more restrictive in terms of allowed activities’"”° 
and have experienced limited increases in human pressures‘. However, there 
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is also evidence to suggest that Indigenous and other mixed-use PAs can be 
more effective because they may have greater local support'**°. 

2. Larger PAs, which have had lower increases in human pressures’, have lower 
deforestation rates as illegal logging deep inside reserves can be logistically 
challenging*'. However, this type of logging is often selective, and thus 
hard to quantify with remotely sensed imagery. Alternatively, larger PAs could 
have higher deforestation rates due to potentially having lower 
budgets per area. 

3. GDP per capita is negatively associated with PA deforestation rates 
because wealthier countries may be better able to fund monitoring and 
enforcement efforts to ensure PA policies are followed. Alternatively, defor- 
estation rates could initially rise with GDP per capita (as a consequence of 
increased resource extraction), and then decline (with increasing funding 
available for conservation)—the environmental Kuznets curve hypothesis”. 
This latter hypothesis would be consistent with the effect of GDP per capita 
on deforestation rates being positive in poorer regions and negative in 
wealthier regions. 

4. Time since establishment (that is, PA age) is negatively associated with defor- 
estation rates because management infrastructure may (rapidly) increase in 
the years after establishment”. 


This list of hypotheses is not exhaustive and, in many cases, we have 
intentionally formed multiple working hypotheses™. Both travel time to nearest 
densely populated area and population density are important measures of potential 
human impacts. However, we did not form a priori hypotheses for these variables 
because it is unclear whether they would have differential impacts on deforestation 
rates in PAs compared with similar unprotected control areas. Similarly, as an 
exploratory analysis, we also tested whether the numbers of threatened and 
non-threatened species were associated with deforestation rates within PAs 
(Supplementary Methods). 


Country indices. We used the data on forest loss in and around PAs to explore 
how area protected and PA effectiveness relate to forest obligate species richness, 
forest carbon stocks and deforestation rates within countries. For each country, we 
calculated: (1) forest obligate vertebrate species richness; (2) forest carbon stocks; 
(3) average forest loss rate; (4) overall PA effectiveness (with respect to limiting 
deforestation); and (5) percentage of forested area protected (Supplementary 
Methods). PA effectiveness at the country level was estimated using the ratio of 
the mean deforestation rate in control areas to the mean deforestation rate of PAs 
within the country. For each country, we refer to this ratio as its PA effectiveness 
score. To simplify interpretation, we restricted our scope to countries with at least 
15 PAs in our main analysis, 5 forest obligate vertebrate species and 10,000 km? of 
forested land. 

We then computed the species threat index, which we defined as the ratio 
of forest species richness to proportion of forested area protected multiplied 
by the PA effectiveness score. Countries with high species threat index scores 
have little protection relative to the number of forest species present. Thus, this 
analysis identifies countries where improvements in PA quantity or quality (that 
is, effectiveness at reducing deforestation rates) may have the most beneficial total 
impact on forest biodiversity. We also calculated variants of the species threat index 
score using the overall forest loss rate and (log-transformed) aboveground forest 
carbon biomass in place of species richness. 

In summary, we computed the following three threat indices, I, for each 
country: 


Species Loss log, )(Carbon) 














T species 34 Loss 34Carbon = 


Prot. x Score Prot. x Score Prot. x Score 


where ‘Species’ is the total number of forest vertebrate species in the country, ‘Prot? 
is the proportion of forested land protected, ‘Score’ (a measure of effectiveness) 

is the total forest loss in control areas divided by total forest loss within PAs, 

‘Loss’ is the estimated annual deforestation rate and ‘Carbon is the aboveground 
forest biomass in units of Gt C. Each index is strictly positive, with higher values 
indicating potential conservation issues because they suggest that the level of 
protection in a country is not commensurate with its biodiversity, deforestation 
pressure, or forest carbon stocks. 


Reporting Summary. Further information on research design is available in the 
Nature Research Reporting Summary linked to this article. 


Data availability 
All data used are publicly available. Sources for the data are given in the Methods 
section. 


Code availability 
Analysis code is available at https://github.com/wolfch2/PA_matching. 
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Extended Data Fig. 1 | The distribution of IUCN categories for the 18,171 PAs in our primary spatial analysis. Protected area categories are: la - ‘Strict 
Nature Reserve; Ib - ‘Wilderness Area, II - ‘National Park,’ III - ‘Natural Monument or Feature,’ IV - ‘Habitat/Species Management Area,’ V - ‘Protected 


Landscape/ Seascape, VI - ‘Protected area with sustainable use of natural resources.’ The protected areas were split into ‘Strict’ (categories I-IV), 
‘Nonstrict’ (categories V-VI), and ‘Unknown’ (any other category). 
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Extended Data Fig. 2 | Net annual forest loss rate within protected areas and in matched control areas. In contrast to the forest loss results, net loss is 
not a true percentage since loss and gain are binary while cover is continuous (see SI Methods for details). Results are grouped by geographic region and 
PA category (IUCN category I-IV: ‘Strict, V-VI: ‘Nonstrict’). Points correspond to median (across PAs) percentage forest loss. Error bar end points are the 
1st and 3rd quartiles for this variable. Forest loss within protected areas has generally been less than in nearby unprotected areas. 
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Extended Data Fig. 3 | Change in the annual forest loss rate associated with the creation of PAs. The change variable is the deforestation rate after minus 
before creation of a PA. Results are grouped by geographic region and PA category (IUCN category I-IV: ‘Strict,’ V-VI: ‘Nonstrict’). Points correspond to 
means, and error bars show standard errors. 
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Extended Data Fig. 4 | Predictors of deforestation rates within protected areas. Each row shows a different predictor variable, and the columns show 


coefficient estimates, standard errors, and FDR-adjusted p 
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geographically. Travel time to nearest densely-populated area was also included as a predictor, but it was found to be non-significant 


spatial variability. Only coefficients with associated p-value less than 0.05 are mapped. 
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Extended Data Fig. 5 | Threatened and non-threatened forest vertebrate species richness. We considered these spatial variables as predictors of 
deforestation within protected areas to explore relationships between PA effectiveness (with respect to limiting deforestation) and biodiversity. 
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Extended Data Fig. 6 | Sensitivity analysis exploring the effect of stricter matching criteria. Medians (center points) and 1* and 3 quartiles (ranges) are 
shown. The first row is for our primary matching dataset (see Fig. 3) based on five classes per continuous matching covariate while the second row shows 
results based on 10 classes per covariate (only 9 were used for travel time - see Supplementary Methods). Overall, the use of stricter matching criteria did 
not appear to considerably alter our results. 
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Extended Data Fig. 7 | Predictors of deforestation rates within protected areas for dataset using stricter matching criteria. Travel time to nearest 


densely 


0.20) was not spatially varying and is not shown in order to parallel our main results (Extended Data Fig. 4). Additionally, 
and strict protection were all found to be constant spatially for this restricted dataset. Only coefficients with associated p-value 


populated area (p 
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PA age 
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less than 0.05 are mapped. 
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Reporting Summary 





Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency 
in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist. 








Statistics 





For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section. 


n/a | Confirmed 





| The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement 








| A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly 








| The statistical test(s) used AND whether they are one- or two-sided 
Only common tests should be described solely by name; describe more complex techniques in the Methods section. 








[| A description of all covariates tested 





A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons 








| A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) 
AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) 








For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted 
Give P values as exact values whenever suitable. 





For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings 








For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes 

















| Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated 








Our web collection on statistics for biologists contains articles on many of the points above. 


Software and code 





Policy information about availability of computer code 





Data collection | No software was used (all data are publicly available). 


Data analysis We carried out the GIS analysis using Google Earth Engine to download most datasets, R v4.0.3 and Python v3.7.3 with GDAL v2.4.0 for 
general raster processing, Julia v1.4.2 for coarsened exact matching, and R v4.0.3 for statistical modeling and data visualization (with 
‘ggplot2’). 


For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and 
reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information. 





Data 





Policy information about availability of data 
All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: 


- Accession codes, unique identifiers, or web links for publicly available datasets 
- A list of figures that have associated raw data 
- A description of any restrictions on data availability 





All data used are publicly available. Sources for the data are given in the methods section. 
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Field-specific reporting 





Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. 


[_] Life sciences [_] Behavioural & social sciences Ecological, evolutionary & environmental sciences 


For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf 





Ecological, evolutionary & environmental sciences study design 





All studies must disclose on these points even when the disclosure is negative. 


Study description This is a global analysis of deforestation rates in protected areas (PAs). The general framework used was a covariate-based matching 
method where deforestation rates in PAs are compared to deforestation rates in matched control areas that serve as a baseline. 


Research sample Not applicable (we used a database of more than 18,000 protected areas). 
Sampling strategy Not applicable (no sampling was conducted). 
Data collection The primary data source is an existing, public database of protected areas (WDPA). Other data sources are also publicly available. 


Timing and spatial scale We considered protected areas (PAs) established in 2000 or earlier for our main analysis and PAs established between 2002 and 
2017 for a before-after comparison. Forest change data range from 2001 to 2018. The PA data are global but restricted to forest 
biomes. 


Data exclusions We did not collect new data for this project. As noted in the SI, we "excluded PAs from our primary analysis that: [1] had point 
(centroid) information only, [2] were exclusively marine, [3] were established after 2000 since the forest loss data range from 2001 to 
2018, [4] had area less than 1 km2 since forest change in small PAs can be hard to estimate accurately, [5] were entirely outside of 
the forest change maps’ common extent, [6] had no land with forest change data within their boundaries, [7] were entirely outside of 
forest biome(s), [8] had less than 30% forest cover across their entire extents, or [9] could not be matched with appropriate control 
areas." 


Reproducibility Not applicable (no experiments were conducted). 


Randomization Not applicable (we used all suitable data in a database of protected areas). 





Blinding Not applicable (this is a global analysis of protected areas). 


Did the study involve field work? C] Yes No 


Reporting for specific materials, systems and methods 





We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, 
system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. 



































Materials & experimental systems Methods 
n/a | Involved in the study n/a | Involved in the study 
l Antibodies ChIP-seq 
l Eukaryotic cell lines Flow cytometry 
m Palaeontology and archaeology |_| MRI-based neuroimaging 

















L Animals and other organisms 





E Human research participants 





l Clinical data 




















l Dual use research of concern 





Content courtesy of Springer Nature, terms of use apply. Rights reserved 


=) 
je) 
e 
= 
= 
(ap) 
= 
(ap) 
VN 
D 
o% 
= 
Q 
=F 
= 
a) 
O 
(@) 
oa 
= 
Ke) 
Za) 
= 
3 
o% 
5 
< 





Terms and Conditions 


Springer Nature journal content, brought to you courtesy of Springer Nature Customer Service Center GmbH (“Springer Nature”). 

Springer Nature supports a reasonable amount of sharing of research papers by authors, subscribers and authorised users (“Users”), for small- 
scale personal, non-commercial use provided that all copyright, trade and service marks and other proprietary notices are maintained. By 
accessing, sharing, receiving or otherwise using the Springer Nature journal content you agree to these terms of use (“Terms”). For these 
purposes, Springer Nature considers academic use (by researchers and students) to be non-commercial. 

These Terms are supplementary and will apply in addition to any applicable website terms and conditions, a relevant site licence or a personal 
subscription. These Terms will prevail over any conflict or ambiguity with regards to the relevant terms, a site licence or a personal subscription 
(to the extent of the conflict or ambiguity only). For Creative Commons-licensed articles, the terms of the Creative Commons license used will 
apply. 

We collect and use personal data to provide access to the Springer Nature journal content. We may also use these personal data internally within 
ResearchGate and Springer Nature and as agreed share it, in an anonymised way, for purposes of tracking, analysis and reporting. We will not 
otherwise disclose your personal data outside the ResearchGate or the Springer Nature group of companies unless we have your permission as 
detailed in the Privacy Policy. 

While Users may use the Springer Nature journal content for small scale, personal non-commercial use, it is important to note that Users may 
not: 


1. use such content for the purpose of providing other users with access on a regular or large scale basis or as a means to circumvent access 
control; 

2. use such content where to do so would be considered a criminal or statutory offence in any jurisdiction, or gives rise to civil liability, or is 
otherwise unlawful; 

3. falsely or misleadingly imply or suggest endorsement, approval , sponsorship, or association unless explicitly agreed to by Springer Nature in 
writing; 

4. use bots or other automated methods to access the content or redirect messages 

5. override any security feature or exclusionary protocol; or 

6. share the content in order to create substitute for Springer Nature products or services or a systematic database of Springer Nature journal 
content. 


In line with the restriction against commercial use, Springer Nature does not permit the creation of a product or service that creates revenue, 
royalties, rent or income from our content or its inclusion as part of a paid for service or for other commercial gain. Springer Nature journal 
content cannot be used for inter-library loans and librarians may not upload Springer Nature journal content on a large scale into their, or any 
other, institutional repository. 

These terms of use are reviewed regularly and may be amended at any time. Springer Nature is not obligated to publish any information or 
content on this website and may remove it or features or functionality at our sole discretion, at any time with or without notice. Springer Nature 
may revoke this licence to you at any time and remove access to any copies of the Springer Nature journal content which have been saved. 

To the fullest extent permitted by law, Springer Nature makes no warranties, representations or guarantees to Users, either express or implied 
with respect to the Springer nature journal content and all parties disclaim and waive any implied warranties or warranties imposed by law, 
including merchantability or fitness for any particular purpose. 

Please note that these rights do not automatically extend to content, data or other material published by Springer Nature that may be licensed 
from third parties. 

If you would like to use or distribute our Springer Nature journal content to a wider audience or on a regular basis or in any other manner not 
expressly permitted by these Terms, please contact Springer Nature at 


onlineservice@springernature.com 





